Oklahoma Tornadoes

library(sf)           

Linking to GEOS 3.11.2, GDAL 3.6.2, PROJ 9.2.0; sf_use_s2() is TRUE

library(rgdal)        

Please note that rgdal will be retired during October 2023,
plan transition to sf/stars/terra functions using GDAL and PROJ
at your earliest convenience.
See https://r-spatial.org/r/2023/05/15/evolution4.html and https://github.com/r-spatial/evolution
rgdal: version: 1.6-7, (SVN revision 1203)
Geospatial Data Abstraction Library extensions to R successfully loaded
Loaded GDAL runtime: GDAL 3.6.2, released 2023/01/02
Path to GDAL shared files: C:/Users/jcolpitt/AppData/Local/R/win-library/4.3/rgdal/gdal
 GDAL does not use iconv for recoding strings.
GDAL binary built with GEOS: TRUE 
Loaded PROJ runtime: Rel. 9.2.0, March 1st, 2023, [PJ_VERSION: 920]
Path to PROJ shared files: C:\Program Files\PostgreSQL\15\share\contrib\postgis-3.3\proj
PROJ CDN enabled: FALSE
Linking to sp version:1.6-1
To mute warnings of possible GDAL/OSR exportToProj4() degradation,
use options("rgdal_show_exportToProj4_warnings"="none") before loading sp or rgdal.

library(ggplot2)      
library(dplyr)        

Attaching package: ‘dplyr’

The following objects are masked from ‘package:stats’:

    filter, lag

The following objects are masked from ‘package:base’:

    intersect, setdiff, setequal, union

library(tidyr)        
library(scales)       
library(RColorBrewer) 
library(units)

udunits database from C:/Users/jcolpitt/AppData/Local/R/win-library/4.3/units/share/udunits/udunits2.xml

library(cowplot)
library(here)

here() starts at C:/Users/jcolpitt/SpatialDataScience/Rprojects/R

okcounty <- st_read(here("data", "ok_counties.shp"), quiet = TRUE)
tpoint <- st_read(here("data", "ok_tornado_point.shp"), quiet = TRUE)
tpath <- st_read(here("data", "ok_tornado_path.shp"), quiet = TRUE)

class(okcounty)

[1] "sf"         "data.frame"

glimpse(okcounty)

Rows: 77
Columns: 8
$ STATEFP  <chr> "40", "40", "40", "40", "40", "40", "40", "40", "40", "40", "40", "40", "40", "40", "40", "40", "40", "40", "40", "40", "40", "40", "40", "40", "40", "40", "40", "40", "40", "40", "40", …
$ COUNTYFP <chr> "077", "025", "011", "107", "105", "153", "001", "053", "059", "073", "129", "015", "097", "033", "147", "041", "109", "057", "133", "093", "031", "063", "075", "099", "141", "083", "069…
$ COUNTYNS <chr> "01101826", "01101800", "01101793", "01101841", "01101840", "01101864", "01101788", "01101814", "01101817", "01101824", "01101852", "01101795", "01101833", "01101804", "01101861", "01101…
$ AFFGEOID <chr> "0500000US40077", "0500000US40025", "0500000US40011", "0500000US40107", "0500000US40105", "0500000US40153", "0500000US40001", "0500000US40053", "0500000US40059", "0500000US40073", "05000…
$ GEOID    <chr> "40077", "40025", "40011", "40107", "40105", "40153", "40001", "40053", "40059", "40073", "40129", "40015", "40097", "40033", "40147", "40041", "40109", "40057", "40133", "40093", "40031…
$ NAME     <chr> "Latimer", "Cimarron", "Blaine", "Okfuskee", "Nowata", "Woodward", "Adair", "Grant", "Harper", "Kingfisher", "Roger Mills", "Caddo", "Mayes", "Cotton", "Washington", "Delaware", "Oklahom…
$ LSAD     <chr> "06", "06", "06", "06", "06", "06", "06", "06", "06", "06", "06", "06", "06", "06", "06", "06", "06", "06", "06", "06", "06", "06", "06", "06", "06", "06", "06", "06", "06", "06", "06", …
$ geometry <POLYGON [°]> POLYGON ((-95.50766 35.0292..., POLYGON ((-103.0025 36.6751..., POLYGON ((-98.6369 36.16489..., POLYGON ((-96.62486 35.4627..., POLYGON ((-95.80982 36.9419..., POLYGON ((-99.6055…

okcounty_sp <- as(okcounty, 'Spatial')
class(okcounty_sp)

[1] "SpatialPolygonsDataFrame"
attr(,"package")
[1] "sp"

okcounty_sf <- st_as_sf(okcounty_sp)
class(okcounty_sf)

[1] "sf"         "data.frame"

ggplot(data = okcounty) +
  geom_sf(fill = NA)

names(tpoint)

 [1] "om"       "yr"       "mo"       "dy"       "date"     "time"     "tz"       "st"       "stf"      "stn"     
[11] "mag"      "inj"      "fat"      "loss"     "closs"    "slat"     "slon"     "elat"     "elon"     "len"     
[21] "wid"      "fc"       "geometry"

tpoint_16_21 <- tpoint %>%
  filter(yr >= 2016 & yr <= 2021) %>%
  select(om, yr, date)

tpath_16_21 <- tpath %>%
  filter(yr >= 2016 & yr <= 2021) %>%
  select(om, yr, date)

ggplot() +
  geom_sf(data = okcounty, fill = NA) +
  geom_sf(data = tpoint_16_21) +
  theme_bw()

ggplot() +
  geom_sf(data = okcounty, fill = NA) +
  geom_sf(data = tpath_16_21, color = "red") +
  theme_void()

ggplot() +
  geom_sf(data = tpoint_16_21,
          aes(color = as.factor(yr))) +
  geom_sf(data = okcounty, fill = NA) +
  scale_color_discrete(name = "Year") +
  coord_sf(datum = NA) +
  theme_void()

ggplot() +
  geom_sf(data = okcounty, 
          fill = NA,
          color = "gray") +
  geom_sf(data = tpoint_16_21,size = 0.75) +
  facet_wrap(vars(yr), ncol = 2) +
  coord_sf(datum = NA) +
  theme_void()

countypnt <- st_join(tpoint_16_21, okcounty)
glimpse(countypnt)

Rows: 434
Columns: 11
$ om       <dbl> 613662, 613675, 613676, 613677, 613678, 613727, 613728, 613729, 613730, 613731, 613763, 613764, 61376…
$ yr       <dbl> 2016, 2016, 2016, 2016, 2016, 2016, 2016, 2016, 2016, 2016, 2016, 2016, 2016, 2016, 2016, 2016, 2016,…
$ date     <chr> "2016-03-23", "2016-03-30", "2016-03-30", "2016-03-30", "2016-03-30", "2016-04-15", "2016-04-15", "20…
$ STATEFP  <chr> "40", "40", "40", "40", "40", "40", "40", "40", "40", "40", "40", "40", "40", "40", "40", "40", "40",…
$ COUNTYFP <chr> "001", "113", "105", "131", "035", "139", "139", "139", "139", "139", "017", "109", "109", "113", "12…
$ COUNTYNS <chr> "01101788", "01101844", "01101840", "01101853", "01101805", "01101857", "01101857", "01101857", "0110…
$ AFFGEOID <chr> "0500000US40001", "0500000US40113", "0500000US40105", "0500000US40131", "0500000US40035", "0500000US4…
$ GEOID    <chr> "40001", "40113", "40105", "40131", "40035", "40139", "40139", "40139", "40139", "40139", "40017", "4…
$ NAME     <chr> "Adair", "Osage", "Nowata", "Rogers", "Craig", "Texas", "Texas", "Texas", "Texas", "Texas", "Canadian…
$ LSAD     <chr> "06", "06", "06", "06", "06", "06", "06", "06", "06", "06", "06", "06", "06", "06", "06", "06", "06",…
$ geometry <POINT [°]> POINT (-94.5042 35.6916), POINT (-96.0151 36.2151), POINT (-95.5523 36.7291), POINT (-95.6384 3…

# convert to a non spatial dataframe
countypnt <- st_drop_geometry(countypnt)
countysum <- countypnt %>%
  group_by(GEOID) %>%
  summarize(tcnt = n())
glimpse(countysum)

Rows: 75
Columns: 2
$ GEOID <chr> "40001", "40005", "40007", "40009", "40011", "40013", "40015", "40017", "40019", "40021", "40023", "4002…
$ tcnt  <int> 6, 3, 4, 8, 1, 4, 10, 5, 7, 5, 3, 12, 10, 5, 5, 1, 7, 9, 7, 8, 2, 2, 2, 3, 9, 5, 4, 4, 2, 1, 4, 3, 9, 1,…

countymap <- okcounty %>%
  left_join(countysum, by = "GEOID") %>%
  replace(is.na(.), 0) %>%
  mutate(area = st_area(okcounty),
         tdens = 10^6 * 10^3 * tcnt / area) %>%
  drop_units()
glimpse(countymap)

Rows: 77
Columns: 11
$ STATEFP  <chr> "40", "40", "40", "40", "40", "40", "40", "40", "40", "40", "40", "40", "40", "40", "40", "40", "40",…
$ COUNTYFP <chr> "077", "025", "011", "107", "105", "153", "001", "053", "059", "073", "129", "015", "097", "033", "14…
$ COUNTYNS <chr> "01101826", "01101800", "01101793", "01101841", "01101840", "01101864", "01101788", "01101814", "0110…
$ AFFGEOID <chr> "0500000US40077", "0500000US40025", "0500000US40011", "0500000US40107", "0500000US40105", "0500000US4…
$ GEOID    <chr> "40077", "40025", "40011", "40107", "40105", "40153", "40001", "40053", "40059", "40073", "40129", "4…
$ NAME     <chr> "Latimer", "Cimarron", "Blaine", "Okfuskee", "Nowata", "Woodward", "Adair", "Grant", "Harper", "Kingf…
$ LSAD     <chr> "06", "06", "06", "06", "06", "06", "06", "06", "06", "06", "06", "06", "06", "06", "06", "06", "06",…
$ tcnt     <dbl> 1, 12, 1, 10, 6, 2, 6, 0, 4, 9, 3, 10, 12, 1, 2, 8, 13, 4, 5, 5, 5, 1, 9, 3, 16, 6, 9, 1, 8, 5, 4, 7,…
$ geometry <POLYGON [°]> POLYGON ((-95.50766 35.0292..., POLYGON ((-103.0025 36.6751..., POLYGON ((-98.6369 36.16489..…
$ area     <dbl> 1890663261, 4766283042, 2427121029, 1657249513, 1503893122, 3227039659, 1493619771, 2597526414, 26851…
$ tdens    <dbl> 0.5289149, 2.5176851, 0.4120108, 6.0340944, 3.9896452, 0.6197631, 4.0170866, 0.0000000, 1.4896951, 3.…

st_write(countymap, dsn = "oktornado.gpkg",
         layer = "countysum",
         delete_dsn = TRUE)

st_write(countymap, dsn = "oktornado.geojson",
         layeroptions = "RFC7946 = YES")

# proud to say the highest tornado density is my home town, Tulsa County
ggplot(data = countymap) +
  geom_sf(aes(fill = tdens)) +
  theme_void()

st_geometry_type(okcounty, by_geometry = FALSE)

[1] POLYGON
18 Levels: GEOMETRY POINT LINESTRING POLYGON MULTIPOINT MULTILINESTRING ... TRIANGLE

okcntrd = st_centroid(countymap)

Warning: st_centroid assumes attributes are constant over geometries

st_geometry_type(okcntrd, by_geometry = FALSE)

[1] POINT
18 Levels: GEOMETRY POINT LINESTRING POLYGON MULTIPOINT MULTILINESTRING ... TRIANGLE

ggplot() +
  geom_sf(data = okcntrd, aes(size = tcnt)) +
  geom_sf(data = okcounty, fill = NA) +
  theme_void()

ggplot(data = countymap) +
  geom_sf(aes(fill = tdens)) +
  scale_fill_distiller(name = expression("Tornadoes/1000 km"^2),
                       palette = "YlOrRd",
                       breaks = pretty_breaks(),
                       direction = 1) +
  theme_void() +
  theme(legend.position = "bottom")

numclas <- 4
qbrks <- seq(0, 1, length.out = numclas + 1)
qbrks

[1] 0.00 0.25 0.50 0.75 1.00

countymap <- countymap %>%
  mutate(tdens_c1 = cut(tdens,
                        breaks = quantile(tdens, breaks = qbrks),
                        include.lowest = T))


ggplot(data = countymap) +
  geom_sf(aes(fill = tdens_c1)) +
  scale_fill_brewer(name = expression("Tornadoes/1000 km"^2),
                    palette = "YlOrRd") +
  theme_void() +
  theme(legend.position = "bottom")

maxcnt <- max(okcntrd$tcnt)
brkpts <- c(0, 2, 5, 10, maxcnt)
okcntrd <- okcntrd %>%
  mutate(tcnt_c1 = cut(tcnt,
                       breaks = brkpts,
                       include.lowest = T))

ggplot(data = okcntrd) +
  geom_sf(aes(size = tcnt_c1)) +
  scale_size_discrete(name="Tornadoes") +
  geom_sf(data = okcounty, fill = NA) +
  theme_void() +
  theme(legend.position = "bottom")

Warning: Using size for a discrete variable is not advised.

ggsave("tornado.png", 
       width = 6, 
       height = 4, 
       dpi = 300)
ggsave("tornado2.png", 
       width = 15, 
       height = 10, 
       units = "cm", 
       dpi = 100)

choropleth <- ggplot(data = countymap) +
  geom_sf(aes(fill = tdens_c1)) +
  scale_fill_brewer(name="Density",   
                    palette = "YlOrRd") +
  theme_void()

gradsymbol <- ggplot(data = okcntrd) +
  geom_sf(aes(size = tcnt_c1)) +
  scale_size_discrete(name="Count") +
  geom_sf(data = okcounty, fill = NA) +
  theme_void()

Warning: Using size for a discrete variable is not advised.

ggsave("choropleth.tiff", 
       plot = choropleth,
       width = 6, 
       height = 4, 
       dpi = 300, 
       compression = "lzw")

ggsave("gradsymbol.tiff",
       plot = gradsymbol,
       width = 6, 
       height = 4, 
       dpi = 300, 
       compression = "lzw")

plot_grid(choropleth, gradsymbol, 
          labels = c("A) Choropleth Map", 
                     "B) Graduated Symbol Map",
                     label_size = 12),
          ncol = 1, 
          hjust = 0, 
          label_x = 0)

tpath_years <- ggplot() +
                  geom_sf(data = tpath_16_21,
                          aes(color = as.factor(yr))) +
                  geom_sf(data = okcounty, fill = NA) +
                  scale_color_discrete(name = "Year") +
                  coord_sf(datum = NA) +
                  theme_void()

tpoint_years <- ggplot() +
                  geom_sf(data = tpoint_16_21,
                          aes(color = as.factor(yr))) +
                  geom_sf(data = okcounty, fill = NA) +
                  scale_color_discrete(name = "Year") +
                  coord_sf(datum = NA) +
                  theme_void()

ggsave("tpath_years.tiff", 
       plot = choropleth,
       width = 6, 
       height = 6, 
       dpi = 300, 
       compression = "lzw")

ggsave("tpoint_years.tiff", 
       plot = choropleth,
       width = 6, 
       height = 6, 
       dpi = 300, 
       compression = "lzw")

plot_grid(tpath_years, tpoint_years, 
          labels = c("A) Tornado Paths", 
                     "B) Tornado Points",
                     label_size = 12),
          ncol = 1, 
          rel_heights = c(1, 1),
          align = "hv",
          vjust = 5)

ggplot(data = countymap) +
  geom_sf(aes(fill = tdens)) +
  scale_fill_distiller(name = expression("Tornadoes/1000 km"^2),
                       palette = "YlOrRd",
                       breaks = pretty_breaks(),
                       direction = 1) +

  theme_void() +
  theme(legend.position = "bottom")

df_join <- countymap %>%
  left_join(countypnt, by = "GEOID") %>%
  replace(is.na(.), 0) %>%
  filter(yr != 0)  

ggplot(data = df_join) +
  geom_sf(aes(fill = tdens)) +
  geom_sf(data = okcounty, 
          fill = NA,
          color = "gray") +
  scale_fill_distiller(name = expression("Tornadoes/1000 km"^2),
                       palette = "YlOrRd",
                       breaks = pretty_breaks(),
                       direction = 1) +
  facet_wrap(vars(yr), ncol = 2) +
  coord_sf(datum = NA) +
  theme_void() +
  theme(legend.position = "bottom")

url <- "https://earthquake.usgs.gov/earthquakes/feed/v1.0/summary/2.5_month.geojson"
earthquakes <- readOGR(url)

Warning: OGR support is provided by the sf and terra packages among othersWarning: OGR support is provided by the sf and terra packages among othersWarning: OGR support is provided by the sf and terra packages among othersWarning: OGR support is provided by the sf and terra packages among othersWarning: OGR support is provided by the sf and terra packages among othersWarning: OGR support is provided by the sf and terra packages among othersWarning: OGR support is provided by the sf and terra packages among others

OGR data source with driver: GeoJSON 
Source: "https://earthquake.usgs.gov/earthquakes/feed/v1.0/summary/2.5_month.geojson", layer: "2.5_month"
with 1947 features
It has 27 fields
Integer64 fields read as strings:  time updated 

eqsf <- st_as_sf(earthquakes)
faults_url <- "https://services.arcgis.com/jIL9msH9OI208GCb/ArcGIS/rest/services/Active_Faults/FeatureServer/0/query?where=1%3D1&objectIds=&time=&geometry=&geometryType=esriGeometryEnvelope&inSR=&spatialRel=esriSpatialRelIntersects&resultType=none&distance=0.0&units=esriSRUnit_Meter&relationParam=&returnGeodetic=false&outFields=*&returnGeometry=true&featureEncoding=esriDefault&multipatchOption=xyFootprint&maxAllowableOffset=&geometryPrecision=&outSR=&defaultSR=&datumTransformation=&applyVCSProjection=false&returnIdsOnly=false&returnUniqueIdsOnly=false&returnCountOnly=false&returnExtentOnly=false&returnQueryGeometry=false&returnDistinctValues=false&cacheHint=false&orderByFields=&groupByFieldsForStatistics=&outStatistics=&having=&resultOffset=&resultRecordCount=&returnZ=false&returnM=false&returnExceededLimitFeatures=true&quantizationParameters=&sqlFormat=none&f=pgeojson&token="
faults <-readOGR(faults_url)

Warning: OGR support is provided by the sf and terra packages among othersWarning: OGR support is provided by the sf and terra packages among othersWarning: OGR support is provided by the sf and terra packages among othersWarning: OGR support is provided by the sf and terra packages among othersWarning: OGR support is provided by the sf and terra packages among othersWarning: OGR support is provided by the sf and terra packages among othersWarning: OGR support is provided by the sf and terra packages among others

OGR data source with driver: GeoJSON 
Source: "https://services.arcgis.com/jIL9msH9OI208GCb/ArcGIS/rest/services/Active_Faults/FeatureServer/0/query?where=1%3D1&objectIds=&time=&geometry=&geometryType=esriGeometryEnvelope&inSR=&spatialRel=esriSpatialRelIntersects&resultType=none&distance=0.0&units=esriSRUnit_Meter&relationParam=&returnGeodetic=false&outFields=*&returnGeometry=true&featureEncoding=esriDefault&multipatchOption=xyFootprint&maxAllowableOffset=&geometryPrecision=&outSR=&defaultSR=&datumTransformation=&applyVCSProjection=false&returnIdsOnly=false&returnUniqueIdsOnly=false&returnCountOnly=false&returnExtentOnly=false&returnQueryGeometry=false&returnDistinctValues=false&cacheHint=false&orderByFields=&groupByFieldsForStatistics=&outStatistics=&having=&resultOffset=&resultRecordCount=&returnZ=false&returnM=false&returnExceededLimitFeatures=true&quantizationParameters=&sqlFormat=none&f=pgeojson&token=", layer: "OGRGeoJSON"
with 16121 features
It has 6 fields

faults_sf <- st_as_sf(faults)

pal <- colorBin(
  palette = "Spectral",
  domain = eqsf$mag,
  reverse = TRUE,
  bins = 5
)

leaflet() %>%
  addTiles() %>%
  setView(-117.841293, 46.195042, 3) %>%
  addProviderTiles(providers$CartoDB, group = "Grayscale") %>% 
  addProviderTiles(providers$Esri.WorldTerrain, group = "Terrain") %>%
  addPolylines(data = faults_sf, group = "vectorData") %>%
  addCircleMarkers(data = eqsf,
                   fillColor = ~pal(mag),
                   radius = ~eqsf$mag * 2,
                   stroke = FALSE,
                   color = "Spectral",
                   fillOpacity = 0.6,
                   popup = paste0(
                     "<strong>Title:</strong> ", eqsf$title,
                     "<br><strong>Magnitude:</strong> ", eqsf$mag,
                     "<br><strong>Intensity:</strong> ", eqsf$mmi,
                     "<br><strong>Sig:</strong> ", eqsf$sig
                   ),
                   group = "vectorData"
  ) %>% 
  addLegend(pal = pal, values = eqsf$mag, position = "bottomleft", title = "Magnitude") %>%
  addLayersControl(overlayGroups = c("vectorData"), baseGroups = c("Terrain", "Grayscale"))